Trabajo preliminar
Inicialización de las redes
Cálculo de las distancias
- Matrices de distancias
- Union de las matrices con la base de viajes
Estadísticas
Cartografía en Español
Cartography in English
Correlaciones
Palma ratio
Valores promediados de los indicadores

Trabajo preliminar

Paquetes R

library(sp)
library(sf)
library(sfnetworks)
library(dodgr)
library(r5r)
library(tidygraph)
library(data.table)
library(dplyr)
library(purrr)
library(TSP)
library(openxlsx)
library(spatstat)
library(raster)
library(stars)
library(terra)
library(ggplot2)
library(ggspatial)
library(knitr)
library(tidyr)
library(stringr)
library(RColorBrewer)
library(foreach)
library(doFuture)
library(RColorBrewer)

Cargando los datos del cuestionario

ZAT <- st_read(dsn = "231031_ZAT_2023", layer = "231031_ZAT_2023") %>% st_transform(4326) %>% filter(UTAM != "N/A")

## Reading layer `231031_ZAT_2023' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota 2023\231031_ZAT_2023' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1215 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY, XYZ
## Bounding box:  xmin: 909732.5 ymin: 904238.5 xmax: 1113716 ymax: 1137247
## z_range:       zmin: 0 zmax: 0
## Projected CRS: MAGNA-SIRGAS / Colombia Bogota zone

ZATCentroids <- st_centroid(ZAT)
 
setwd("C:/Users/hugot/Documents/Rennes 2/Calcul distance SIG/Reprise perso/RMD Bogota")

UTAM <- st_read(dsn = "Data", layer = "EMU2019_area")

## Reading layer `EMU2019_area' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 132 features and 57 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 569984 ymin: 492795.1 xmax: 626083.9 ymax: 557354.6
## Projected CRS: WGS 84 / UTM zone 18N

EF <- read.xlsx("Data/facteurs d'émissions.xlsx")
limites_bogota <- st_read(dsn = "Data", layer = "Localidad_Municipio_2017") %>% # Municipalities boundaries
  st_transform(4326)

## Reading layer `Localidad_Municipio_2017' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 564193.1 ymin: 483995.6 xmax: 634113.2 ymax: 570281.8
## Projected CRS: WGS 1984 UTM, Zone 18 North, Meter

setwd("C:/Users/hugot/Documents/Uniandes/EODH2023/EODH 2023 oficial/Publicacion1804/EODH/02_Base datos procesada/XLSX")

viajes <- read.xlsx("d. Modulo viajes.xlsx")
Hogares <- read.xlsx("a. Modulo hogares.xlsx")
Personas <- read.xlsx("c. Modulo personas.xlsx")
#Etapas <- read.xlsx("e. Modulo etapas.xlsx")
#Vehiculos <- read.xlsx("b. Modulo vehiculos.xlsx")

setwd("C:/Users/hugot/Documents/Uniandes/EODH2023")

#Esta base es necesaria para tener la ubicación del hogar
Hogares_loc <- read.xlsx(xlsxFile = "ANEXO 2_BASE DE DATOS EODH CON FACTORES DE EXPANSIÓN.xlsx", sheet = "hogar A-C-X-CF")
#Estas base son necesarias para deducir si el lugar de origen y destino de cada viaje es hogar u otro
Personas_viaj_no_procesada <- read.xlsx(xlsxFile = "ANEXO 2_BASE DE DATOS EODH CON FACTORES DE EXPANSIÓN.xlsx", sheet = "modulo D-E")
Viajes_no_procesada <- read.xlsx(xlsxFile = "ANEXO 2_BASE DE DATOS EODH CON FACTORES DE EXPANSIÓN.xlsx", sheet = "viajes")

Cargando los datos SHP

setwd("C:/Users/hugot/Documents/Rennes 2/Calcul distance SIG/Reprise perso/RMD Bogota")

Routes <- st_read(dsn = "Data", layer = "Reseau_routier")

## Reading layer `Reseau_routier' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 100036 features and 12 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 565362.9 ymin: 486923.8 xmax: 633899.8 ymax: 569868
## Projected CRS: WGS 84 / UTM zone 18N

Routes <- st_cast(Routes, "LINESTRING") %>% st_transform(4326)
Transmi <- st_read(dsn = "Data", layer = "Transmilenio") %>% st_transform(4326)

## Reading layer `Transmilenio' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 13 features and 7 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 588057.1 ymin: 503694.6 xmax: 606100.3 ymax: 527398.5
## Projected CRS: WGS 84 / UTM zone 18N

Transmi_stops <- st_read(dsn = "Data", layer = "estaciones-de-transmilenio") %>% st_transform(4326)

## Reading layer `estaciones-de-transmilenio' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 140 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.20572 ymin: 4.55681 xmax: -74.04358 ymax: 4.768743
## Geodetic CRS:  WGS 84

Procesos preliminares

Se definen las categorías de modos que se usarán dentro del modelo de distancias

modos <- data.frame(modo_principal_desagrupado = unique(viajes$modo_principal_desagrupado),
                    modo_principal_typo = c("Auto", "Bicicleta", "A pie", "A pie", "Intermunicipal", "TransMilenio", "Transporte publico individual", "SITP Zonal", "Moto", "Alimentador", "Moto", "Auto", "Transporte publico individual", "Auto", "Transporte publico individual", "Bicitaxi", "SITP Zonal", "Transporte Escolar", "Auto", "Transporte informal", "Otro", "Otro", "Otro", "Transporte informal", "Auto", "Otro", "Bicicleta", "Cable", "Otro", "Patineta", "Bicicleta", "Otro", "Otro"),
                    modo_principal_comparado = c("Auto", "Bicicleta", "A pie", "A pie", "Transporte publico", "Transporte publico", "Taxi / Carro por aplicación", "Transporte publico", "Moto", "Transporte publico", "Moto", "Auto", "Taxi / Carro por aplicación", "Auto", "Taxi / Carro por aplicación", "Transporte informal", "Transporte publico", "Transporte Escolar", "Auto", "Transporte informal", "Otro", "Otro", "Otro", "Transporte informal", "Auto", "Otro", "Bicicleta", "Transporte publico", "Otro", "Otro", "Bicicleta", "Otro", "Otro"),
                    modo_principal_LATS = c("Private car", "Bike", "Walking", "Walking", "Regular bus", "BRT", "Taxi", "Regular bus", "Moto", "Regular bus", "Moto", "Private car", "Taxi", "Private car", "Taxi", "Paratransit", "BRT", "School bus", "Private car", "Paratransit", "Other", "Other", "Other", "Paratransit", "Private car", "Other", "Bike", "BRT", "Other", "Other", "Bike", "Other", "Other"))

viajes <- viajes %>% left_join(modos, by = "modo_principal_desagrupado")

Excluimos los viajes entrando o saliendo del area conformado por el Distrito de Bogotá y los 20 municipios que fueron encuestados. Para ello, se conservan solo los viajes que empiezan y terminan en las ZAT para los cuales tenemos un archivo SHP con su geometría y a los cuales les corresponde una UTAM. A la fecha, se usa el archivo con las 1,215 ZAT de 2023, y se reduce el número de viajes en la base de 100,174 a 96,054 (-4%).

#UTAM_presentes <- data.frame(UTAM = unique(c(Hogares$cod_utam_hg, viajes$utam_des, viajes$utam_ori))) %>% filter(UTAM != "No aplica")
#ZAT_presentes <- data.frame(ZAT = unique(c(Hogares$zat_hg, viajes$zat_des, viajes$zat_ori)))
#write.xlsx(UTAM_presentes, "UTAM_EODH_2023.xlsx")
#write.xlsx(ZAT_presentes, "ZAT_EODH_2023.xlsx")

viajes <- viajes[viajes$zat_des %in% ZAT$ZAT & viajes$zat_ori %in% ZAT$ZAT,]

Necesitamos recuperar las coordenadas de los hogares desde la versión de trabajo de la encuesta.

Hogares <- Hogares %>% 
  left_join(Hogares_loc[,c("i7a-Latitude", "i7a-Longitude", "KEY")], by = c("key_hg" = "KEY")) %>%
  st_as_sf(coords = c("i7a-Longitude", "i7a-Latitude"), crs = 4326) 

Hogares_loc <- Hogares_loc[,c("i7a-Latitude", "i7a-Longitude", "cod_dane","KEY")] %>%
  st_as_sf(coords = c("i7a-Longitude", "i7a-Latitude"), crs = 4326) 

#st_write(Hogares_loc, "Hogares_loc_mza.shp")

#Definición de la ventana
zoom_level <- 9.5
lon_span <- 360 / 2^zoom_level
lat_span <- 360 / 2^zoom_level

zoom_to_Bogota <- c(-74.12, 4.72)  # Bogota
lon_bounds_Bogota <- c(zoom_to_Bogota[1] - lon_span / 2, zoom_to_Bogota[1] + lon_span / 2)
lat_bounds_Bogota <- c(zoom_to_Bogota[2] - lat_span / 2, zoom_to_Bogota[2] + lat_span / 2)

UTAM <- UTAM %>% st_transform(4326)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, fill = "white")+
  geom_sf(data = Hogares_loc, col = "red", size = 0.5)+
  labs(title = "Ubicación de los hogares de la muestra", 
       subtitle = "Bogotá DC + 20 municipios", 
       caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, Secretaría Distrital de Movilidad, 2023", 
       col = "UPZ") +
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(x = "", y = "") +
  theme(legend.position = "right",
        legend.text = element_text(size=10),
        legend.title = element_text(size=15),
        plot.title = element_text(size=15),
        plot.caption = element_text(size = 7, face = "italic", hjust = 0, vjust = 12),
        legend.text.align = 1)+
  annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(0.8, "cm"))

Finalmente necesitamos identificar para cada viaje si sus puntos de origen y destino son el hogar de la persona u otro. Para hacerlo, necesitamos buscar en la versión de trabajo de la base de viajes la variable que contiene esta información para el lugar de destino únicamente (lastimosamente no la hay para el lugar de origen). Luego, para el lugar de origen del primer viaje del día, buscamos esta información en la version de trabajo de la base de viajes. Este procedimiento busca aumentar la precisión de los viajes que comienzan o terminan en el hogar

#Lugar destino
viajes <- viajes %>% left_join(Viajes_no_procesada[,c("d24","KEY")], by = c("key_viaje" = "KEY"))
colnames(viajes)[colnames(viajes) == "d24"] <- "lugar_destino"
viajes$lugar_destino[viajes$lugar_destino == 1] <- "Hogar"
viajes$lugar_destino[viajes$lugar_destino == 2] <- "Otro lugar"

#Lugar origen (funciona solo para el primer viaje del día)
Personas_viaj_no_procesada$KEY_2 <- gsub(pattern = "person_d_e", replacement = "B/person_b", Personas_viaj_no_procesada$KEY)
Personas <- Personas %>% left_join(Personas_viaj_no_procesada[,c("d7", "KEY_2")], by = c("key_persona" = "KEY_2"))
colnames(Personas)[colnames(Personas) == "d7"] <- "lugar_ini_dia"
Personas$lugar_ini_dia[Personas$lugar_ini_dia == 1] <- "Hogar"
viajes <- viajes %>% left_join(Personas[,c("cod_per", "lugar_ini_dia")], by = c("cod_pers" = "cod_per"))
viajes$lugar_origen <- "Otro lugar"
viajes$lugar_origen[viajes$lugar_ini_dia == "Hogar" & viajes$orden_vj == 1] <- "Hogar"

#Control (no debe haber viajes Hogar --> Hogar)
viajes$test_origen_destino[viajes$lugar_origen == "Hogar" & viajes$lugar_destino == "Hogar"] <- "HOGAR_HOGAR"
viajes$test_origen_destino[viajes$lugar_origen == "Hogar" & viajes$lugar_destino == "Otro lugar"] <- "HOGAR_Otro"
viajes$test_origen_destino[viajes$lugar_origen == "Otro lugar" & viajes$lugar_destino == "Hogar"] <- "Otro_HOGAR"
viajes$test_origen_destino[viajes$lugar_origen == "Otro lugar" & viajes$lugar_destino == "Otro lugar"] <- "Otro_Otro"
print(table(viajes$test_origen_destino))

## 
## HOGAR_HOGAR  HOGAR_Otro  Otro_HOGAR   Otro_Otro 
##           1       45499       45721        4833

#Encontramos un (1) viaje "Hogar --> Hogar". Quitamos todos los viajes de la persona que lo realiza del registro de viajes
excl <- viajes$cod_pers[viajes$test_origen_destino == "HOGAR_HOGAR"]
viajes <- viajes[!(viajes$cod_pers %in% excl),]

Inicialización de las redes

Red vial

# Ponderación de los segmentos y restricción al principal componente conexo
net_mot <- weight_streetnet(Routes, wt_profile = "motorcar", type_col = "type")
net_mot <- net_mot[net_mot$component == 1,]
net_cycl <- weight_streetnet(Routes, wt_profile = "bicycle", type_col = "type")
net_cycl <- net_cycl[net_cycl$component == 1,]
net_walk <- weight_streetnet(Routes, wt_profile = "foot", type_col = "type")
net_walk <- net_walk[net_walk$component == 1,]

Transmilenio

Transmi <- st_cast(Transmi, "LINESTRING")
Transmi_stops <- Transmi_stops %>% st_transform(4326)

# Inicialización
net_transmi = as_sfnetwork(Transmi, directed = FALSE) %>%
  st_transform(4326) %>%
  activate("edges") %>%
  mutate(weight = edge_length())

# Subdivisión de ejes
net_transmi = convert(net_transmi, to_spatial_subdivision)

# Inclusión (blend) de las estaciones a la red
net_transmi = st_network_blend(net_transmi, Transmi_stops)

# Cálculo de la nueva longitud de cada eje
net_transmi <- net_transmi %>% activate("edges") %>% mutate(weight = edge_length())

SITP (R5R)

#Allocating 8 GB of memory to Java (2 GB is a minimum, adjust according to the hardware)
options(java.parameters = "-Xmx12G")
#gc()

#Initiating Java
#rJava::.jinit()
#rJava::.jcall("java.lang.System", "S", "getProperty", "java.version")

#Setting up the multimodal network. Requires a .zip archive containing the GTFS and a .pbf file with the road network.
#It is quite long the first time (~15 min), then it reloads the previously created network.
r5r_core <- setup_r5(data_path = "C:/Users/hugot/Documents/Rennes 2/Calcul distance SIG/Reprise perso/RMD Bogota/Data")

# set inputs
mode <- c("WALK", "TRANSIT")
departure_datetime <- as.POSIXct("13-05-2019 14:00:00",
                                 format = "%d-%m-%Y %H:%M:%S")

#Putting coordinates in the appropriate format for r5r
HogaresR5R <- cbind(as.data.frame(st_coordinates(Hogares)), Hogares$cod_hog)
colnames(HogaresR5R) <- c("lon", "lat", "id")

ZATR5R <- cbind(as.data.frame(st_coordinates(ZATCentroids)), ZATCentroids$ZAT)
colnames(ZATR5R) <- c("lon", "lat", "id")

Cálculo de las distancias

Matrices de distancias

Se calcula la longitud de cada viaje de acuerdo a los siguiente:

Se le asigna un modo principal
Se escoge una red de acuerdo al modo principal (ver tabla adjunta) para la asignación del viaje origen-destino.
Se calculan las rutas desde y hasta los centroides de las ZAT de origen y destino (si éste es cualquiera (caso general), o desde/hasta el hogar cuando aplica.

Modo principal	Red empleada
Transmilenio	Transmilenio
SITP Zonal	SITP
Alimentador	SITP
A pie	Red vial (ponderación peatón)
Bicicleta	Red vial (ponderación ciclista)
Cualquier otro	Red vial (ponderación automovilista)

print(Sys.time())

#5 min

#LINEA RECTA (para sustituir NA)
# 1 --> 2 (Hogar --> ZAT) 
mat_straight_12 <- as.data.frame(st_distance(x = st_geometry(Hogares), y = st_geometry(ZATCentroids)))
rownames(mat_straight_12) <- Hogares$cod_hog
colnames(mat_straight_12) <- ZAT$ZAT
mat_straight_12_table <- as.data.frame(as.table(as.matrix(mat_straight_12)))
colnames(mat_straight_12_table) <- c("cod_hog", "zat_des", "dist_straight")

# 2 --> 1 (ZAT --> Hogar) 
mat_straight_21 <- as.data.frame(st_distance(x = st_geometry(ZATCentroids), y = st_geometry(Hogares)))
rownames(mat_straight_21) <- ZAT$ZAT
colnames(mat_straight_21) <- Hogares$cod_hog
mat_straight_21_table <- as.data.frame(as.table(as.matrix(mat_straight_21)))
colnames(mat_straight_21_table) <- c("zat_ori", "cod_hog", "dist_straight")

# 2 --> 2 (ZAT --> ZAT) 
mat_straight_22 <- as.data.frame(st_distance(x = st_geometry(ZATCentroids), y = st_geometry(ZATCentroids)))
rownames(mat_straight_22) <- ZAT$ZAT
colnames(mat_straight_22) <- ZAT$ZAT
mat_straight_22_table <- as.data.frame(as.table(as.matrix(mat_straight_22)))
colnames(mat_straight_22_table) <- c("zat_ori", "zat_des", "dist_straight")

#TRANSMILENIO
# 1 --> 2 (Hogar --> ZAT) 
mat_transmi_12 <- st_network_cost(net_transmi, from = st_geometry(Hogares), to = st_geometry(ZATCentroids))
rownames(mat_transmi_12) <- Hogares$cod_hog
colnames(mat_transmi_12) <- ZAT$ZAT
mat_transmi_12_table <- as.data.frame(as.table(mat_transmi_12))
colnames(mat_transmi_12_table) <- c("cod_hog", "zat_des", "dist_transmi")

# 2 --> 1 (ZAT --> Hogar) 
mat_transmi_21 <- st_network_cost(net_transmi, from = st_geometry(ZATCentroids), to = st_geometry(Hogares))
rownames(mat_transmi_21) <- ZAT$ZAT
colnames(mat_transmi_21) <- Hogares$cod_hog
mat_transmi_21_table <- as.data.frame(as.table(mat_transmi_21))
colnames(mat_transmi_21_table) <- c("zat_ori", "cod_hog", "dist_transmi")

# 2 --> 2 (ZAT --> ZAT) 
mat_transmi_22 <- st_network_cost(net_transmi, from = st_geometry(ZATCentroids), to = st_geometry(ZATCentroids))
rownames(mat_transmi_22) <- ZAT$ZAT
colnames(mat_transmi_22) <- ZAT$ZAT
mat_transmi_22_table <- as.data.frame(as.table(mat_transmi_22))
colnames(mat_transmi_22_table) <- c("zat_ori", "zat_des", "dist_transmi")

#MOTORIZADO
# 1 --> 2 (Hogar --> ZAT) (10 secs)
mat_mot_12 <- dodgr_dists(graph = net_mot, from = as.data.frame(st_coordinates(Hogares)), to = as.data.frame(st_coordinates(ZATCentroids)), shortest = TRUE)
rownames(mat_mot_12) <- Hogares$cod_hog
colnames(mat_mot_12) <- ZAT$ZAT
mat_mot_12[is.na(mat_mot_12)] <- mat_straight_12[is.na(mat_mot_12)]
mat_mot_12_table <- as.data.frame(as.table(mat_mot_12))
colnames(mat_mot_12_table) <- c("cod_hog", "zat_des", "dist_car")

# 2 --> 1 (ZAT --> Hogar) (10 secs)
mat_mot_21 <- dodgr_dists(graph = net_mot, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(Hogares)), shortest = TRUE)
rownames(mat_mot_21) <- ZAT$ZAT
colnames(mat_mot_21) <- Hogares$cod_hog
mat_mot_21[is.na(mat_mot_21)] <- mat_straight_21[is.na(mat_mot_21)]
mat_mot_21_table <- as.data.frame(as.table(mat_mot_21))
colnames(mat_mot_21_table) <- c("zat_ori", "cod_hog", "dist_car")

# 2 --> 2 (ZAT --> ZAT) (10 secs)
mat_mot_22 <- dodgr_dists(graph = net_mot, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(ZATCentroids)), shortest = TRUE)
rownames(mat_mot_22) <- ZAT$ZAT
colnames(mat_mot_22) <- ZAT$ZAT
mat_mot_22[is.na(mat_mot_22)] <- mat_straight_22[is.na(mat_mot_22)]
mat_mot_22_table <- as.data.frame(as.table(mat_mot_22))
colnames(mat_mot_22_table) <- c("zat_ori", "zat_des", "dist_car")

#BICI
# 1 --> 2 (Hogar --> ZAT) (10 secs)
mat_cycl_12 <- dodgr_dists(graph = net_cycl, from = as.data.frame(st_coordinates(Hogares)), to = as.data.frame(st_coordinates(ZATCentroids)), shortest = TRUE)
rownames(mat_cycl_12) <- Hogares$cod_hog
colnames(mat_cycl_12) <- ZAT$ZAT
mat_cycl_12[is.na(mat_cycl_12)] <- mat_straight_12[is.na(mat_cycl_12)]
mat_cycl_12_table <- as.data.frame(as.table(mat_cycl_12))
colnames(mat_cycl_12_table) <- c("cod_hog", "zat_des", "dist_bike")

# 2 --> 1 (ZAT --> Hogar) (10 secs)
mat_cycl_21 <- dodgr_dists(graph = net_cycl, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(Hogares)), shortest = TRUE)
rownames(mat_cycl_21) <- ZAT$ZAT
colnames(mat_cycl_21) <- Hogares$cod_hog
mat_cycl_21[is.na(mat_cycl_21)] <- mat_straight_21[is.na(mat_cycl_21)]
mat_cycl_21_table <- as.data.frame(as.table(mat_cycl_21))
colnames(mat_cycl_21_table) <- c("zat_ori", "cod_hog", "dist_bike")

# 2 --> 2 (ZAT --> ZAT) (10 secs)
mat_cycl_22 <- dodgr_dists(graph = net_cycl, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(ZATCentroids)), shortest = TRUE)
rownames(mat_cycl_22) <- ZAT$ZAT
colnames(mat_cycl_22) <- ZAT$ZAT
mat_cycl_22[is.na(mat_cycl_22)] <- mat_straight_22[is.na(mat_cycl_22)]
mat_cycl_22_table <- as.data.frame(as.table(mat_cycl_22))
colnames(mat_cycl_22_table) <- c("zat_ori", "zat_des", "dist_bike")

#WALK
# 1 --> 2 (Hogar --> ZAT) (10 secs)
mat_walk_12 <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(Hogares)), to = as.data.frame(st_coordinates(ZATCentroids)), shortest = TRUE)
rownames(mat_walk_12) <- Hogares$cod_hog
colnames(mat_walk_12) <- ZAT$ZAT
mat_walk_12[is.na(mat_walk_12)] <- mat_straight_12[is.na(mat_walk_12)]
mat_walk_12_table <- as.data.frame(as.table(mat_walk_12))
colnames(mat_walk_12_table) <- c("cod_hog", "zat_des", "dist_walk")

# 2 --> 1 (ZAT --> Hogar) (10 secs)
mat_walk_21 <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(Hogares)), shortest = TRUE)
rownames(mat_walk_21) <- ZAT$ZAT
colnames(mat_walk_21) <- Hogares$cod_hog
mat_walk_21[is.na(mat_walk_21)] <- mat_straight_21[is.na(mat_walk_21)]
mat_walk_21_table <- as.data.frame(as.table(mat_walk_21))
colnames(mat_walk_21_table) <- c("zat_ori", "cod_hog", "dist_walk")

# 2 --> 2 (ZAT --> ZAT) (10 secs)
mat_walk_22 <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(ZATCentroids)), shortest = TRUE)
rownames(mat_walk_22) <- ZAT$ZAT
colnames(mat_walk_22) <- ZAT$ZAT
mat_walk_22[is.na(mat_walk_22)] <- mat_straight_22[is.na(mat_walk_22)]
mat_walk_22_table <- as.data.frame(as.table(mat_walk_22))
colnames(mat_walk_22_table) <- c("zat_ori", "zat_des", "dist_walk")

print(Sys.time())

Union de las matrices con la base de viajes

colnames(viajes)[colnames(viajes) == "cod_hg"] <- "cod_hog"

#Separando la base en 3 sub-bases de acuerdo al lugar de origen y destino

viajes_12 <- viajes %>% filter(test_origen_destino == "HOGAR_Otro")
viajes_21 <- viajes %>% filter(test_origen_destino == "Otro_HOGAR")
viajes_22 <- viajes %>% filter(test_origen_destino == "Otro_Otro")

#Unión por columnas múltiples (15 min)

print(Sys.time())

viajes_12 <- merge(viajes_12, mat_mot_12_table, by = c("cod_hog", "zat_des")) 
viajes_21 <- merge(viajes_21, mat_mot_21_table, by = c("cod_hog", "zat_ori")) 
viajes_22 <- merge(viajes_22, mat_mot_22_table, by = c("zat_ori", "zat_des")) 

viajes_12 <- merge(viajes_12, mat_cycl_12_table, by = c("cod_hog", "zat_des")) 
viajes_21 <- merge(viajes_21, mat_cycl_21_table, by = c("cod_hog", "zat_ori")) 
viajes_22 <- merge(viajes_22, mat_cycl_22_table, by = c("zat_ori", "zat_des")) 

viajes_12 <- merge(viajes_12, mat_walk_12_table, by = c("cod_hog", "zat_des")) 
viajes_21 <- merge(viajes_21, mat_walk_21_table, by = c("cod_hog", "zat_ori")) 
viajes_22 <- merge(viajes_22, mat_walk_22_table, by = c("zat_ori", "zat_des")) 

viajes_12 <- merge(viajes_12, mat_transmi_12_table, by = c("cod_hog", "zat_des")) 
viajes_21 <- merge(viajes_21, mat_transmi_21_table, by = c("cod_hog", "zat_ori")) 
viajes_22 <- merge(viajes_22, mat_transmi_22_table, by = c("zat_ori", "zat_des")) 

viajes_12 <- merge(viajes_12, mat_straight_12_table, by = c("cod_hog", "zat_des")) 
viajes_21 <- merge(viajes_21, mat_straight_21_table, by = c("cod_hog", "zat_ori")) 
viajes_22 <- merge(viajes_22, mat_straight_22_table, by = c("zat_ori", "zat_des"))

print(Sys.time())

## 27 HORAS

print(Sys.time())

#1-> 2 (Hogar -> ZAT)
viajes_12_sitp <- viajes_12[viajes_12$modo_principal_typo == "SITP Zonal",]
viajes_12_origen <- cbind(id = viajes_12_sitp$cod_hog, HogaresR5R[match(viajes_12_sitp$cod_hog, HogaresR5R$id), 1:2, drop = FALSE])
viajes_12_destino <- cbind(id = viajes_12_sitp$zat_des, ZATR5R[match(viajes_12_sitp$zat_des, ZATR5R$id), 1:2, drop = TRUE])

viajes_12_sitp$dist_sitp = 0
for(i in 1:nrow(viajes_12_sitp)){
  #print(i)
  det <- detailed_itineraries(r5r_core = r5r_core,
                                origins = viajes_12_origen[i,],
                                destinations = viajes_12_destino[i,],
                                mode = mode,
                                departure_datetime = departure_datetime,
                                max_walk_time = 180,
                                max_trip_duration = 300,
                                shortest_path = TRUE)
  if(!is.na(det$total_distance[1])){
      det <- det %>% filter(mode == "BUS")
      viajes_12_sitp$dist_sitp[i] <- sum(det$distance)
  }
}

print(Sys.time())
write.xlsx(viajes_12_sitp, "viajes_12_sitp.xlsx")

#2-> 1 (ZAT -> Hogar)
viajes_21_sitp <- viajes_21[viajes_21$modo_principal_typo == "SITP Zonal",]
viajes_21_origen <- cbind(id = viajes_21_sitp$zat_ori, ZATR5R[match(viajes_21_sitp$zat_ori, ZATR5R$id), 1:2, drop = TRUE])
viajes_21_destino <-  cbind(id = viajes_21_sitp$cod_hog, HogaresR5R[match(viajes_21_sitp$cod_hog, HogaresR5R$id), 1:2, drop = FALSE])

viajes_21_sitp$dist_sitp = 0
for(i in 1:nrow(viajes_21_sitp)){
  #print(i)
  det <- detailed_itineraries(r5r_core = r5r_core,
                                origins = viajes_21_origen[i,],
                                destinations = viajes_21_destino[i,],
                                mode = mode,
                                departure_datetime = departure_datetime,
                                max_walk_time = 180,
                                max_trip_duration = 300,
                                shortest_path = TRUE)
  if(!is.na(det$total_distance[1])){
      det <- det %>% filter(mode == "BUS")
      viajes_21_sitp$dist_sitp[i] <- sum(det$distance)
  }
}

print(Sys.time())
write.xlsx(viajes_21_sitp, "viajes_21_sitp.xlsx")

#2-> 2 (ZAT -> ZAT)
viajes_22_sitp <- viajes_22[viajes_22$modo_principal_typo == "SITP Zonal",]
viajes_22_origen <- cbind(id = viajes_22_sitp$zat_ori, ZATR5R[match(viajes_22_sitp$zat_ori, ZATR5R$id), 1:2, drop = TRUE])
viajes_22_destino <-  cbind(id = viajes_22_sitp$cod_hog, HogaresR5R[match(viajes_22_sitp$cod_hog, HogaresR5R$id), 1:2, drop = FALSE])

viajes_22_sitp$dist_sitp = 0
for(i in 1:nrow(viajes_22_sitp)){
  #print(i)
  det <- detailed_itineraries(r5r_core = r5r_core,
                                origins = viajes_22_origen[i,],
                                destinations = viajes_22_destino[i,],
                                mode = mode,
                                departure_datetime = departure_datetime,
                                max_walk_time = 180,
                                max_trip_duration = 300,
                                shortest_path = TRUE)
  if(!is.na(det$total_distance[1])){
      det <- det %>% filter(mode == "BUS")
      viajes_22_sitp$dist_sitp[i] <- sum(det$distance)
  }
}

print(Sys.time())
write.xlsx(viajes_22_sitp, "viajes_22_sitp.xlsx")

viajes_12_sitp <- read.xlsx("viajes_12_sitp.xlsx")
viajes_21_sitp <- read.xlsx("viajes_21_sitp.xlsx")
viajes_22_sitp <- read.xlsx("viajes_22_sitp.xlsx")


viajes_12 <- viajes_12 %>%
  left_join(viajes_12_sitp[,c("cod_vj", "dist_sitp")], by = c("cod_vj"))
viajes_21 <- viajes_21 %>%
  left_join(viajes_21_sitp[,c("cod_vj", "dist_sitp")], by = c("cod_vj"))
viajes_22 <- viajes_22 %>%
  left_join(viajes_22_sitp[,c("cod_vj", "dist_sitp")], by = c("cod_vj"))

viajes_dist <- rbind(viajes_12, viajes_21, viajes_22)

viajes_dist$Distancia <- viajes_dist$dist_car

viajes_dist$Distancia[viajes_dist$modo_principal_typo == "TransMilenio"] <- viajes_dist$dist_transmi[viajes_dist$modo_principal_typo == "TransMilenio"]

viajes_dist$Distancia[viajes_dist$modo_principal_typo == "SITP Zonal"] <- viajes_dist$dist_sitp[viajes_dist$modo_principal_typo == "SITP Zonal"]

viajes_dist$Distancia[viajes_dist$modo_principal_typo == "A pie"] <- viajes_dist$dist_walk[viajes_dist$modo_principal_typo == "A pie"]

viajes_dist$Distancia[viajes_dist$modo_principal_typo %in% c("Bicicleta")] <- viajes_dist$dist_bike[viajes_dist$modo_principal_typo %in% c("Bicicleta")]

#viajes_dist$Distancia <- as.numeric(viajes_dist$Distancia)

viajes_dist <- read.xlsx("viajes_dist_no_corr.xlsx")

Para los viajes en Transmilenio, dado que la distancia desde el punto de partida hasta la estación, y desde la estación hasta el punto de llegada, pueden ser largas, calculamos además la distancia por la red vial con ponderación peatón entre estos puntos, y la vamos sumando a la distancia calculada por la red troncal. El cálculo se realiza entonces con dodgr. Lo mismo sucede para los viajes con el SITP Zonal.

plot(net_transmi)

#Conversion de nudos de la red en formato sf. Incluye 140 estaciones + 25 nudos correspondientes a extremidades de ejes que no pudieron ser evitados en la construcción de la red
acceso_transmi <- net_transmi %>%
  activate("nodes") %>%
  st_geometry() %>%
  st_as_sf()

#Matriz de distancias a pie entre hogares y estaciones
mat_walk_hog_transmi <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(Hogares)), to = as.data.frame(st_coordinates(acceso_transmi)), shortest = TRUE)
rownames(mat_walk_hog_transmi) <- Hogares$cod_hog
colnames(mat_walk_hog_transmi) <- rownames(acceso_transmi)

#Búsqueda de la estación de Transmilenio más cercana a cada hogar
distancia <- apply(mat_walk_hog_transmi[,1:dim(mat_walk_hog_transmi)[2]], FUN = min, MARGIN = 1)
name <- data.frame(name = colnames(mat_walk_hog_transmi)[apply(mat_walk_hog_transmi, 1, function(x) which(x == min(x))[1])])
estacion_mas_cerca_hogar <- cbind(name, distancia)
estacion_mas_cerca_hogar$cod_hog <- as.numeric(rownames(estacion_mas_cerca_hogar))

#Matriz de distancias a pie entre centroides de ZAT y estaciones
mat_walk_centroid_transmi <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(acceso_transmi)), shortest = TRUE)
rownames(mat_walk_centroid_transmi) <- ZAT$ZAT
colnames(mat_walk_centroid_transmi) <- rownames(acceso_transmi)

#Búsqueda de la estación de Transmilenio más cercana a cada centroide de ZAT
distancia <- apply(mat_walk_centroid_transmi[,1:dim(mat_walk_centroid_transmi)[2]], FUN = min, MARGIN = 1)
name <- data.frame(name = colnames(mat_walk_centroid_transmi)[apply(mat_walk_centroid_transmi, 1, function(x) which(x == min(x))[1])])
estacion_mas_cerca_centroid <- cbind(name, distancia)
estacion_mas_cerca_centroid$ZAT <- as.numeric(rownames(estacion_mas_cerca_centroid))

#Adición de la distancia hasta y desde Transmilenio en los viajes con este modo
Hog_TM <- viajes_dist %>% left_join(estacion_mas_cerca_hogar[,c("cod_hog", "distancia")], by = "cod_hog") %>%
  group_by(cod_hog) %>%
  summarize(dist_Hog_TM = unique(distancia), dist_TM_Hog = unique(distancia))

ZAT_TM <- viajes_dist %>% left_join(estacion_mas_cerca_centroid[,c("ZAT", "distancia")], by = c("zat_ori" = "ZAT")) %>%
  group_by(zat_ori) %>%
  summarize(dist_ZAT_TM = unique(distancia))

TM_ZAT <- viajes_dist %>% left_join(estacion_mas_cerca_centroid[,c("ZAT", "distancia")], by = c("zat_des" = "ZAT")) %>%
  group_by(zat_des) %>%
  summarize(dist_TM_ZAT = unique(distancia))

viajes_dist <- viajes_dist %>% 
  left_join(Hog_TM, by = "cod_hog") %>%
  left_join(ZAT_TM, by = "zat_ori") %>%
  left_join(TM_ZAT, by = "zat_des")

viajes_dist$dist_hasta_TM[viajes_dist$test_origen_destino == "HOGAR_Otro"] <- viajes_dist$dist_Hog_TM[viajes_dist$test_origen_destino == "HOGAR_Otro"]
viajes_dist$dist_hasta_TM[viajes_dist$test_origen_destino %in% c("Otro_HOGAR", "Otro_Otro")] <- viajes_dist$dist_ZAT_TM[viajes_dist$test_origen_destino %in% c("Otro_HOGAR", "Otro_Otro")]
viajes_dist$dist_desde_TM[viajes_dist$test_origen_destino == "Otro_HOGAR"] <- viajes_dist$dist_TM_Hog[viajes_dist$test_origen_destino == "Otro_HOGAR"]
viajes_dist$dist_desde_TM[viajes_dist$test_origen_destino %in% c("HOGAR_Otro", "Otro_Otro")] <- viajes_dist$dist_TM_ZAT[viajes_dist$test_origen_destino %in% c("HOGAR_Otro", "Otro_Otro")]

#Se trata de una solución provisional para calcular la distancia a la estación del SITP más cercana mientras resolvemos el problema con r5r.

#Paraderos oficiales del SITP. Incluye 8 517 paraderos del SITP
paraderos_SITP <- st_read(dsn = ".", layer = "Paraderos_SITP") %>% st_transform(4326)

## Reading layer `Paraderos_SITP' from data source 
##   `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota 2023' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 8517 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 587059 ymin: 486191.3 xmax: 610279.6 ymax: 532670.2
## Projected CRS: WGS 84 / UTM zone 18N

plot(st_geometry(paraderos_SITP))

#Matriz de distancias a pie entre hogares y paraderos
mat_walk_hog_SITP <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(Hogares)), to = as.data.frame(st_coordinates(paraderos_SITP)), shortest = TRUE)
rownames(mat_walk_hog_SITP) <- Hogares$cod_hog
colnames(mat_walk_hog_SITP) <- rownames(paraderos_SITP)

#Búsqueda del paradero de SITP más cercano a cada hogar
distancia <- apply(mat_walk_hog_SITP[,1:dim(mat_walk_hog_SITP)[2]], FUN = min, MARGIN = 1)
name <- data.frame(name = colnames(mat_walk_hog_SITP)[apply(mat_walk_hog_SITP, 1, function(x) which(x == min(x))[1])])
estacion_mas_cerca_hogar <- cbind(name, distancia)
estacion_mas_cerca_hogar$cod_hog <- as.numeric(rownames(estacion_mas_cerca_hogar))

#Matriz de distancias a pie entre centroides de ZAT y paraderos
mat_walk_centroid_SITP <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(paraderos_SITP)), shortest = TRUE)
rownames(mat_walk_centroid_SITP) <- ZAT$ZAT
colnames(mat_walk_centroid_SITP) <- rownames(paraderos_SITP)

#Búsqueda del paradero de SITP más cercano a cada centroide de ZAT
distancia <- apply(mat_walk_centroid_SITP[,1:dim(mat_walk_centroid_SITP)[2]], FUN = min, MARGIN = 1)
name <- data.frame(name = colnames(mat_walk_centroid_SITP)[apply(mat_walk_centroid_SITP, 1, function(x) which(x == min(x))[1])])
estacion_mas_cerca_centroid <- cbind(name, distancia)
estacion_mas_cerca_centroid$ZAT <- as.numeric(rownames(estacion_mas_cerca_centroid))

#Adición de la distancia hasta y desde el SITP en los viajes con este modo
Hog_SITP <- viajes_dist %>% left_join(estacion_mas_cerca_hogar[,c("cod_hog", "distancia")], by = "cod_hog") %>%
  group_by(cod_hog) %>%
  summarize(dist_Hog_SITP = unique(distancia), dist_SITP_Hog = unique(distancia))

ZAT_SITP <- viajes_dist %>% left_join(estacion_mas_cerca_centroid[,c("ZAT", "distancia")], by = c("zat_ori" = "ZAT")) %>%
  group_by(zat_ori) %>%
  summarize(dist_ZAT_SITP = unique(distancia))

SITP_ZAT <- viajes_dist %>% left_join(estacion_mas_cerca_centroid[,c("ZAT", "distancia")], by = c("zat_des" = "ZAT")) %>%
  group_by(zat_des) %>%
  summarize(dist_SITP_ZAT = unique(distancia))

viajes_dist <- viajes_dist %>% 
  left_join(Hog_SITP, by = "cod_hog") %>%
  left_join(ZAT_SITP, by = "zat_ori") %>%
  left_join(SITP_ZAT, by = "zat_des")

viajes_dist$dist_hasta_SITP[viajes_dist$test_origen_destino == "HOGAR_Otro"] <- viajes_dist$dist_Hog_SITP[viajes_dist$test_origen_destino == "HOGAR_Otro"]
viajes_dist$dist_hasta_SITP[viajes_dist$test_origen_destino %in% c("Otro_HOGAR", "Otro_Otro")] <- viajes_dist$dist_ZAT_SITP[viajes_dist$test_origen_destino %in% c("Otro_HOGAR", "Otro_Otro")]
viajes_dist$dist_desde_SITP[viajes_dist$test_origen_destino == "Otro_HOGAR"] <- viajes_dist$dist_SITP_Hog[viajes_dist$test_origen_destino == "Otro_HOGAR"]
viajes_dist$dist_desde_SITP[viajes_dist$test_origen_destino %in% c("HOGAR_Otro", "Otro_Otro")] <- viajes_dist$dist_SITP_ZAT[viajes_dist$test_origen_destino %in% c("HOGAR_Otro", "Otro_Otro")]

Hallamos la distancia promedio desde cada red hasta los puntos de la matriz (con resolución 200m).

comparacion <- data.frame(red = c("Transmilenio - Hogar","Transmilenio - ZAT","SITP - Hogar", "SITP - ZAT"),
                          dist_promedio_desde_matriz = c(mean(Hog_TM$dist_Hog_TM),
                                                         mean(ZAT_TM$dist_ZAT_TM),
                                                         mean(Hog_SITP$dist_Hog_SITP),
                                                         mean(ZAT_SITP$dist_ZAT_SITP)),
                          dist_mediana_desde_matriz = c(median(Hog_TM$dist_Hog_TM),
                                                         median(ZAT_TM$dist_ZAT_TM),
                                                         median(Hog_SITP$dist_Hog_SITP),
                                                         median(ZAT_SITP$dist_ZAT_SITP)))

colnames(comparacion) <- c("Red", "Promedio", "Mediana")

kable(comparacion, caption = "Comparación de las distancias entre los Hogares, los centroides de ZAT y las redes de transporte colectivo", format.args = list(big.mark = ","))

Comparación de las distancias entre los Hogares, los centroides de ZAT y las redes de transporte colectivo
Red	Promedio	Mediana
Transmilenio - Hogar	5,671.487	1,941.5376
Transmilenio - ZAT	5,408.001	1,912.4224
SITP - Hogar	2,917.074	213.0079
SITP - ZAT	2,779.874	249.1494

Para los trayectos intrazonas, dos metodos fueron probados en la base de la EODH 2019:

El método “del CEREMA”. La distancia de cada viaje es la media raiz cuadrada del area de la zona.
El método “de la media”. Se sortean 10 puntos en cada ZAT y calcula la distancia entre cada par de origen-destino por la red vial con ponderación peatón (ya que 810 de 1255 viajes internos (65%) se realizan con este modo). Luego se asigna al viaje la distancia correspondiente al promedio de los valores calculados excluyendo los ceros.

En nuestro caso, escogemos el método “de la media” que tiende a dar estimaciones superiores al método “del CEREMA”.

Para el caso de los viajes intrazonas en Transmilenio, el criterio para detectarlos es que la distancia por la red de Transmilenio sea nula.

#Seleccion de las ZAT que tienen trayectos internos
viajes_intra <- viajes_dist[viajes_dist$zat_ori == viajes_dist$zat_des & viajes_dist$Distancia == 0,c("zat_ori", "zat_des", "cod_vj", "modo_principal_typo")]
viajes_intra$id_row <- 1:nrow(viajes_intra)
ZAT_intra <- ZAT[ZAT$ZAT %in% as.numeric(unique(viajes_intra$zat_des)),]

#Tomando muestra en cada ZAT seleccionada
#Tamaño de la muestra
s <- 10
d <- nrow(ZAT_intra)

set.seed(1)
ZATSampled <- st_sample(ZAT_intra, rep(s,d))

ZATSampled <- st_sf(ZATSampled) %>%
  st_join(ZAT[,"ZAT"],
          join = st_intersects)

ZATSampled$row <- rep(c(1:s),d)

#Calculando la matriz de distancias por la red vial, ponderación peatón
mat_intra <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(ZATSampled)), to = as.data.frame(st_coordinates(ZATSampled)), shortest = TRUE)

rownames(mat_intra) <- ZATSampled$ZAT
colnames(mat_intra) <- ZATSampled$ZAT

comb <- function(z) {
  merge(viajes_intra, z, by = c("zat_ori","zat_des")) %>% 
    arrange(id_row) %>%
    dplyr::select(c(-zat_ori,-zat_des, -cod_vj, -modo_principal_typo, -id_row))
}

doFuture::registerDoFuture()
future::plan(multisession, workers = 8)
options(future.globals.maxSize= 1048576000)

viajes2 <- foreach(u = 1:s, .combine = 'cbind') %:%
  foreach(v = 1:s, .combine = 'cbind') %dofuture% {
    subset_i = seq(from = u, to = s*d, by = s)
    subset_j = seq(from = v, to = s*d, by = s)
    z <- mat_intra[subset_i,subset_j]
    colnames(z) <- ZAT_intra$ZAT
    rownames(z) <- ZAT_intra$ZAT
    z <- data.frame(rownames(as.matrix(z))[row(as.matrix(z))], 
               colnames(as.matrix(z))[col(as.matrix(z))], 
               c(as.matrix(z)))
    colnames(z) <- c("zat_ori", "zat_des", paste0("dist",(u-1)*s+(v-1)))
    x <- comb(z)
  }

viajes_intra <- cbind(viajes_intra,viajes2)

#índices de las columnas con distancia interna igual a cero
index_zero <- rep("VIDE",s)
for(i in 1:s){
    index_zero[i] <- paste0("dist",(i-1)*s+(i-1))
}

#Distancia interna = promedio de las distancias excluyendo los ceros.
viajes_intra_non_zero <- viajes_intra[,!(colnames(viajes_intra) %in% index_zero)]
viajes_intra$dist_mean_non_zeros <- rowMeans(viajes_intra_non_zero[,6:95])

#Sustituyendo los ceros en la base de distancias
viajes_dist$Distancia[viajes_dist$zat_ori == viajes_dist$zat_des & viajes_dist$Distancia == 0] <- viajes_intra$dist_mean_non_zeros

Identificamos los viajes de ida y vuelta. Un viaje de ida y vuelta se identifica porque cumple las siguientes características:

Se buscan tales viajes dentro del conjunto de viajes realizados por una misma persona identificada por su cod_hog.
Para cada persona, creamos una submuestra correspondiente a los potenciales viajes “de ida”, que son los que originan en su casa, y otra submuesta correspondiente a los potenciales viajes “de vuelta”.
Se unen las dos submuestras por cod_hog y se comparan el modo de los dos viajes (deben coincidir para que sea ida y vuelta), el lugar de destino del viaje de ida y el lugar de origen del viaje de regreso (deben coincidir para que sea ida y vuelta).
Finalmente, se crea una variable dummy en la base de viajes completa indicando si cada viaje es parte de una ida y vuelta o no.

Nota: dado que no pudimos identificar todos los viajes que inician en el hogar, sino solo el primero del día, este algoritmo solo puede detectar máximo una ida y vuelta por persona

#Esta etapa es necesaria en el marco del control de velocidad de la siguiente etapa. Se define como ida y vuelta un trayecto casa -> cierto lugar -> casa realizado con el mismo modo de transporte. El lugar de destino del primer viaje y el lugar de origen del segundo deben coincidir. Si el modo de transporte cambia, no lo definimos como ida y vuelta.

ida <- viajes_dist[viajes_dist$lugar_origen == "Hogar", c("cod_pers", "cod_vj", "lugar_destino", "modo_principal_typo")]
vuelta <- viajes_dist[viajes_dist$lugar_destino == "Hogar", c("cod_pers", "cod_vj", "lugar_origen", "modo_principal_typo")]
ida_vuelta <- na.omit(ida %>% left_join(vuelta, by = "cod_pers"))
ida_vuelta$is_ida_vuelta <- (ida_vuelta$lugar_destino == ida_vuelta$lugar_origen) & (ida_vuelta$modo_principal_typo.x == ida_vuelta$modo_principal_typo.y)
ida <- ida_vuelta[,c("cod_vj.x", "is_ida_vuelta")]
colnames(ida)[colnames(ida) == "cod_vj.x"] <- "cod_vj"
ida <- ida %>% group_by(cod_vj) %>% summarize(is_ida_vuelta = as.logical(max(is_ida_vuelta)))
vuelta <- ida_vuelta[,c("cod_vj.y", "is_ida_vuelta")]
colnames(vuelta)[colnames(vuelta) == "cod_vj.y"] <- "cod_vj"
vuelta <- vuelta %>% group_by(cod_vj) %>% summarize(is_ida_vuelta = as.logical(max(is_ida_vuelta)))
viajes_dist <- viajes_dist %>% left_join(ida, by = c("cod_vj"))
viajes_dist <- viajes_dist %>% left_join(vuelta, by = c("cod_vj"))
viajes_dist$is_ida_vuelta.x[is.na(viajes_dist$is_ida_vuelta.x)] <- FALSE
viajes_dist$is_ida_vuelta.y[is.na(viajes_dist$is_ida_vuelta.y)] <- FALSE
viajes_dist$is_ida_vuelta <- viajes_dist$is_ida_vuelta.x | viajes_dist$is_ida_vuelta.y
viajes_dist$is_ida_vuelta.x <- NULL
viajes_dist$is_ida_vuelta.y <- NULL

#velocidad en km/h
viajes_dist$velocidad <- 60*viajes_dist$Distancia/(1000*viajes_dist$duracion_min)

#8139 viajes a pie son demasiado rápidos (23%)
viajes_dist$Distancia[viajes_dist$modo_principal_typo == "A pie" & viajes_dist$velocidad > 10] <- 1000*viajes_dist$duracion_min[viajes_dist$modo_principal_typo == "A pie" & viajes_dist$velocidad > 10]*10/60

#250 viajes en bici son demasiado rápidos (5%)
viajes_dist$Distancia[viajes_dist$modo_principal_typo == "Bicicleta" & viajes_dist$velocidad > 30] <- 1000*viajes_dist$duracion_min[viajes_dist$modo_principal_typo == "Bicicleta" & viajes_dist$velocidad > 30]*30/60

#963 viajes en modos motorizados son demasiado rápidos (2%)
viajes_dist$Distancia[!(viajes_dist$modo_principal_typo %in% c("A pie","Bicicleta")) & viajes_dist$velocidad > 50] <- 1000*viajes_dist$duracion_min[!(viajes_dist$modo_principal_typo %in% c("A pie","Bicicleta")) & viajes_dist$velocidad > 50]*50/60

#se aplica una ultima correccion: los viajes de ida y vuelta deben tener la misma distancia. Se corrige la distancia de los viajes demasiado rápidos, pero en el caso de los viajes de ida y vuelta, dado que esperamos que la distancia coincida en la ida y el regreso, le vamos a asignar el maximo entre los dos valores.
ida_vuelta <- viajes_dist[viajes_dist$is_ida_vuelta == TRUE,] %>% 
  group_by(cod_pers, modo_principal_typo) %>%
  summarise(Distancia_AR = max(Distancia))
viajes_dist <- viajes_dist %>% left_join(ida_vuelta, by = c("cod_pers", "modo_principal_typo"))
viajes_dist$Distancia[viajes_dist$is_ida_vuelta == TRUE] <- viajes_dist$Distancia_AR[viajes_dist$is_ida_vuelta == TRUE]
viajes_dist$Distancia_AR <- NULL
viajes_dist$velocidad <- NULL

viajes_dist$Distancia_incl_caminata <- viajes_dist$Distancia

viajes_dist$Distancia_incl_caminata[viajes_dist$modo_principal_typo == "TransMilenio"] <- viajes_dist$Distancia[viajes_dist$modo_principal_typo == "TransMilenio"] + viajes_dist$dist_hasta_TM[viajes_dist$modo_principal_typo == "TransMilenio"] + viajes_dist$dist_desde_TM[viajes_dist$modo_principal_typo == "TransMilenio"]

viajes_dist$Distancia_incl_caminata[viajes_dist$modo_principal_typo == "SITP Zonal"] <- viajes_dist$Distancia[viajes_dist$modo_principal_typo == "SITP Zonal"] + viajes_dist$dist_hasta_SITP[viajes_dist$modo_principal_typo == "SITP Zonal"] + viajes_dist$dist_desde_SITP[viajes_dist$modo_principal_typo == "SITP Zonal"]

Por fin podemos crear la base final con las variables correspondientes a distancia, emisiones de gases de efecto invernadero y contaminantes aéreos.

viajes_dist <- viajes_dist %>% left_join(EF, by = c("modo_principal_typo" = "modo_principal"))

viajes_dist$`CO2-eq` <- viajes_dist$`CO2-eq`*viajes_dist$Distancia/1000
viajes_dist$CO <- viajes_dist$CO*viajes_dist$Distancia/1000
viajes_dist$NOx <- viajes_dist$NOx*viajes_dist$Distancia/1000
viajes_dist$SO2 <- viajes_dist$SO2*viajes_dist$Distancia/1000
viajes_dist$COV <- viajes_dist$COV*viajes_dist$Distancia/1000
viajes_dist$PM.2.5 <- viajes_dist$PM.2.5*viajes_dist$Distancia/1000
viajes_dist$PM.10 <- viajes_dist$PM.10*viajes_dist$Distancia/1000

#write.xlsx(viajes_dist, paste0("viajes_dist_final_con_caminata_hasta_TM.xlsx"), overwrite = TRUE)

Estadísticas

viajes_dist <- read.xlsx("viajes_dist_final_con_caminata_hasta_TM.xlsx")

estadisticas <- viajes_dist %>% 
  group_by(modo_principal_typo) %>% 
  summarize(numero_viajes = round(sum(fexp_vj),0), 
            mpkt = round(sum(fexp_vj*Distancia/1000000000),1),
            mpkt_incl_caminata = round(sum(fexp_vj*Distancia_incl_caminata/1000000000),1),
            reparto_modal = round(100*sum(fexp_vj)/sum(viajes_dist$fexp_vj),1), 
            reparto_modal_pkt = round(100*sum(fexp_vj*Distancia_incl_caminata)/sum(viajes_dist$fexp_vj*viajes_dist$Distancia_incl_caminata),1), 
            `CO2-eq`= round(sum(fexp_vj*`CO2-eq`/1000000),1),
            dist_por_viaj = round(sum(fexp_vj*Distancia_incl_caminata/1000)/sum(fexp_vj),1))
  
estadisticas <- estadisticas[order(estadisticas$numero_viajes, decreasing = TRUE),]

total <- data.frame("TOTAL", 
                   sum(estadisticas$numero_viajes),
                   sum(estadisticas$mpkt),
                   sum(estadisticas$mpkt_incl_caminata),
                   sum(estadisticas$reparto_modal),
                   sum(estadisticas$reparto_modal_pkt),
                   sum(estadisticas$`CO2-eq`),
                   round(sum(estadisticas$numero_viajes*estadisticas$dist_por_viaj)/sum(estadisticas$numero_viajes),1))

colnames(total) <- colnames(estadisticas)

estadisticas <- rbind(estadisticas, total)

colnames(estadisticas) <- c("Modo principal","Número de viajes","Millón PKT (solo tramo principal)", "Millón PKT (total incluyendo caminata inicial y final)","% Número de viajes","% PKT","Emisiones GEI (tCO2-eq)", "Distancia media por viaje (km)")

kable(estadisticas, caption = "Estadísticas básicas", format.args = list(big.mark = ","))

Estadísticas básicas
Modo principal	Número de viajes	Millón PKT (solo tramo principal)	Millón PKT (total incluyendo caminata inicial y final)	% Número de viajes	% PKT	Emisiones GEI (tCO2-eq)	Distancia media por viaje (km)
A pie	6,134,682	10.8	10.8	37.4	9.1	0.0	1.8
SITP Zonal	2,237,192	24.6	25.9	13.6	22.0	1,469.7	11.6
TransMilenio	2,008,743	26.3	34.0	12.3	28.8	105.8	16.9
Auto	1,916,020	16.0	16.0	11.7	13.6	2,433.2	8.4
Moto	1,053,826	10.9	10.9	6.4	9.3	834.5	10.4
Bicicleta	1,037,486	5.7	5.7	6.3	4.9	0.0	5.5
Transporte publico individual	670,709	4.3	4.3	4.1	3.6	695.8	6.4
Intermunicipal	432,446	4.1	4.1	2.6	3.5	293.7	9.6
Transporte Escolar	293,703	1.9	1.9	1.8	1.6	96.4	6.4
Otro	250,446	2.1	2.1	1.5	1.8	106.0	8.3
Alimentador	178,527	1.1	1.1	1.1	0.9	9.2	6.1
Transporte informal	119,618	0.8	0.8	0.7	0.7	64.0	7.0
Bicitaxi	47,601	0.2	0.2	0.3	0.2	10.5	4.8
Patineta	12,677	0.1	0.1	0.1	0.1	0.2	5.9
Cable	2,269	0.0	0.0	0.0	0.0	0.1	3.5
TOTAL	16,395,945	108.9	117.9	99.9	100.1	6,119.1	7.2

#write.xlsx(estadisticas, "resumen_distancias_emisiones_con_caminata_hasta_TM.xlsx", overwrite = TRUE)

estadisticas <- viajes_dist %>% 
  group_by(modo_principal_comparado) %>% 
  summarize(numero_viajes = round(sum(fexp_vj),0), 
            mpkt = round(sum(fexp_vj*Distancia/1000000000),1),
            mpkt_incl_caminata = round(sum(fexp_vj*Distancia_incl_caminata/1000000000),1),
            reparto_modal = round(100*sum(fexp_vj)/sum(viajes_dist$fexp_vj),1), 
            reparto_modal_pkt = round(100*sum(fexp_vj*Distancia_incl_caminata)/sum(viajes_dist$fexp_vj*viajes_dist$Distancia_incl_caminata),1), 
            `CO2-eq`= round(sum(fexp_vj*`CO2-eq`/1000000),1),
            dist_por_viaj = round(sum(fexp_vj*Distancia_incl_caminata/1000)/sum(fexp_vj),1))
  
estadisticas <- estadisticas[order(estadisticas$numero_viajes, decreasing = TRUE),]

total <- data.frame("TOTAL", 
                   sum(estadisticas$numero_viajes),
                   sum(estadisticas$mpkt),
                   sum(estadisticas$mpkt_incl_caminata),
                   sum(estadisticas$reparto_modal),
                   sum(estadisticas$reparto_modal_pkt),
                   sum(estadisticas$`CO2-eq`),
                   round(sum(estadisticas$numero_viajes*estadisticas$dist_por_viaj)/sum(estadisticas$numero_viajes),1))

colnames(total) <- colnames(estadisticas)

estadisticas <- rbind(estadisticas, total)

colnames(estadisticas) <- c("Modo principal","Número de viajes","Millón PKT (solo tramo principal)", "Millón PKT (total incluyendo caminata inicial y final)","% Número de viajes","% PKT","Emisiones GEI (tCO2-eq)", "Distancia media por viaje (km)")

kable(estadisticas, caption = "Estadísticas básicas", format.args = list(big.mark = ","))

Estadísticas básicas
Modo principal	Número de viajes	Millón PKT (solo tramo principal)	Millón PKT (total incluyendo caminata inicial y final)	% Número de viajes	% PKT	Emisiones GEI (tCO2-eq)	Distancia media por viaje (km)
A pie	6,134,682	10.8	10.8	37.4	9.1	0.0	1.8
Transporte publico	4,859,177	56.1	65.1	29.6	55.2	1,878.5	13.4
Auto	1,916,020	16.0	16.0	11.7	13.6	2,433.2	8.4
Moto	1,053,826	10.9	10.9	6.4	9.3	834.5	10.4
Bicicleta	1,037,486	5.7	5.7	6.3	4.9	0.0	5.5
Taxi / Carro por aplicación	670,709	4.3	4.3	4.1	3.6	695.8	6.4
Transporte Escolar	293,703	1.9	1.9	1.8	1.6	96.4	6.4
Otro	263,123	2.1	2.1	1.6	1.8	106.1	8.2
Transporte informal	167,219	1.1	1.1	1.0	0.9	74.5	6.4
TOTAL	16,395,945	108.9	117.9	99.9	100.0	6,119.0	7.2

#write.xlsx(estadisticas, "comparado_distancias_emisiones_con_caminata_hasta_TM.xlsx", overwrite = TRUE)

estadisticas <- viajes_dist %>% 
  group_by(modo_principal_LATS) %>% 
  summarize(numero_viajes = round(sum(fexp_vj),0), 
            mpkt = round(sum(fexp_vj*Distancia/1000000000),1),
            mpkt_incl_caminata = round(sum(fexp_vj*Distancia_incl_caminata/1000000000),1),
            reparto_modal = round(100*sum(fexp_vj)/sum(viajes_dist$fexp_vj),1), 
            reparto_modal_pkt = round(100*sum(fexp_vj*Distancia_incl_caminata)/sum(viajes_dist$fexp_vj*viajes_dist$Distancia_incl_caminata),1), 
            `CO2-eq`= round(sum(fexp_vj*`CO2-eq`/1000000),1),
            dist_por_viaj = round(sum(fexp_vj*Distancia_incl_caminata/1000)/sum(fexp_vj),1))
  
estadisticas <- estadisticas[order(estadisticas$numero_viajes, decreasing = TRUE),]

total <- data.frame("TOTAL", 
                   sum(estadisticas$numero_viajes),
                   sum(estadisticas$mpkt),
                   sum(estadisticas$mpkt_incl_caminata),
                   sum(estadisticas$reparto_modal),
                   sum(estadisticas$reparto_modal_pkt),
                   sum(estadisticas$`CO2-eq`),
                   round(sum(estadisticas$numero_viajes*estadisticas$dist_por_viaj)/sum(estadisticas$numero_viajes),1))

colnames(total) <- colnames(estadisticas)

estadisticas <- rbind(estadisticas, total)

colnames(estadisticas) <- c("Modo principal","Número de viajes","Millón PKT (solo tramo principal)", "Millón PKT (total incluyendo caminata inicial y final)","% Número de viajes","% PKT","Emisiones GEI (tCO2-eq)", "Distancia media por viaje (km)")

kable(estadisticas, caption = "Estadísticas básicas", format.args = list(big.mark = ","))

Estadísticas básicas
Modo principal	Número de viajes	Millón PKT (solo tramo principal)	Millón PKT (total incluyendo caminata inicial y final)	% Número de viajes	% PKT	Emisiones GEI (tCO2-eq)	Distancia media por viaje (km)
Walking	6,134,682	10.8	10.8	37.4	9.1	0.0	1.8
Regular bus	2,751,356	28.6	29.9	16.8	25.4	1,704.2	10.9
BRT	2,107,820	27.5	35.2	12.9	29.8	174.3	16.7
Private car	1,916,020	16.0	16.0	11.7	13.6	2,433.2	8.4
Moto	1,053,826	10.9	10.9	6.4	9.3	834.5	10.4
Bike	1,037,486	5.7	5.7	6.3	4.9	0.0	5.5
Taxi	670,709	4.3	4.3	4.1	3.6	695.8	6.4
Transporte Escolar	293,703	1.9	1.9	1.8	1.6	96.4	6.4
Other	263,123	2.1	2.1	1.6	1.8	106.1	8.2
Paratransit	167,219	1.1	1.1	1.0	0.9	74.5	6.4
TOTAL	16,395,944	108.9	117.9	100.0	100.0	6,119.0	7.2

#write.xlsx(estadisticas, "LATS_distancias_emisiones_con_caminata_hasta_TM.xlsx", overwrite = TRUE)

Cartografía en Español

Finalmente, producimos mapas por zonas (UTAM)

centroides_limites_bogota <- st_centroid(limites_bogota) %>%
  mutate(lon = st_coordinates(.)[,1],
         lat = st_coordinates(.)[,2])

centroides_limites_bogota$MUNI <- centroides_limites_bogota$NOMBRE
centroides_limites_bogota$MUNI[centroides_limites_bogota$NOMBR_MUNI == "BOGOTA"] <- ""

centroides_limites_bogota$LOCA <- centroides_limites_bogota$NOMBRE
centroides_limites_bogota$LOCA[centroides_limites_bogota$NOMBR_MUNI != "BOGOTA"] <- ""

centroides_limites_bogota$DISPLAY <- centroides_limites_bogota$NOMBRE
centroides_limites_bogota$DISPLAY[!((centroides_limites_bogota$NOMBR_MUNI == "BOGOTA"  & centroides_limites_bogota$NOMBRE %in% c("BOSA", "SANTA FE", "CIUDAD BOLIVAR", "TEUSAQUILLO", "FONTIBON", "ENGATIVA", "SUBA", "USAQUEN", "USME", "SAN CRISTOBAL", "CHAPINERO", "KENNEDY")) | (centroides_limites_bogota$NOMBR_MUNI != "BOGOTA"))] <- ""

centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "SUBA"] <- 4.76
centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "SOACHA"] <- 4.57
centroides_limites_bogota$lon[centroides_limites_bogota$NOMBRE == "CHAPINERO"] <- -74.027
centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "USAQUEN"] <- 4.74
centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "CHAPINERO"] <- 4.68
centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "TEUSAQUILLO"] <- 4.65
centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "KENNEDY"] <- 4.635
centroides_limites_bogota$DISPLAY[centroides_limites_bogota$NOMBRE == "CIUDAD BOLIVAR"] <- "CIUDAD \nBOLIVAR"

ZAT_Poblacion <- Personas %>%
  group_by(zat_hg) %>%
  summarize(pop = sum(`fexp_per>5años`))

ZAT <- ZAT %>% left_join(ZAT_Poblacion, by = c("ZAT" = "zat_hg"))
ZAT$pop[is.na(ZAT$pop)] <- 0

UTAM <- UTAM %>% st_transform(4326)

Dist_ZAT <- viajes_dist[,c("zat_hg", "modo_principal_comparado", "Distancia_incl_caminata", "fexp_vj")]%>%
  left_join(ZAT[,c("ZAT")], by = c("zat_hg" = "ZAT")) %>%
  group_by(zat_hg) %>%
  summarize(dist = sum(fexp_vj*Distancia_incl_caminata/1000))

ZAT <- ZAT %>% left_join(Dist_ZAT, by = c("ZAT" = "zat_hg"))

ZAT$dist[is.na(ZAT$dist)] <- 0

UTAM_indic <- ZAT %>% 
  st_drop_geometry() %>%
  group_by(UTAM) %>% 
  summarise(dist = sum(dist), 
           pop = sum(pop), 
           dist_capita = dist/pop)

UTAM <- UTAM %>% left_join(UTAM_indic, by = "UTAM")

#Colores
bks_total <- round(as.numeric(quantile(UTAM$dist_capita, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- c(bks_total,max(UTAM$dist_capita))

UTAM$binned_dist_capita <- cut(UTAM$dist_capita,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = binned_dist_capita))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Límite de municipio fuera de Bogotá DC \nLímite de localidad dentro de Bogotá DC"), color = "grey")+
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 4.5, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(6,"Purples")) +
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(title = "", fill = "Distancia cotidiana promedio per \ncápita todos modos incluidos (km)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, Secretaría Distrital de Movilidad, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 90),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

Dist_ZAT <- viajes_dist[viajes_dist$modo_principal_comparado == "A pie",c("zat_hg", "modo_principal_comparado", "Distancia_incl_caminata", "fexp_vj")]%>%
  left_join(ZAT[,c("ZAT")], by = c("zat_hg" = "ZAT")) %>%
  group_by(zat_hg) %>%
  summarize(dist_walk = sum(fexp_vj*Distancia_incl_caminata/1000))

ZAT <- ZAT %>% left_join(Dist_ZAT, by = c("ZAT" = "zat_hg"))

ZAT$dist_walk[is.na(ZAT$dist_walk)] <- 0

UTAM_indic <- ZAT %>% 
  st_drop_geometry() %>%
  group_by(UTAM) %>% 
  summarise(dist_walk = sum(dist_walk), 
           walk_intensity = sum(dist_walk)/sum(pop))

UTAM <- UTAM %>% left_join(UTAM_indic, by = "UTAM")

#Colores
bks_walk <- round(as.numeric(quantile(UTAM$walk_intensity, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- c(bks_walk,max(UTAM$walk_intensity))

UTAM$binned_walk_intensity <- cut(UTAM$walk_intensity,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = binned_walk_intensity))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Límite de municipio fuera de Bogotá DC \nLímite de localidad dentro de Bogotá DC"), color = "grey")+
  #geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 4.5, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(6,"Reds")) +
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(title = "", fill = "Distancia cotidiana promedio \nper cápita caminando (km)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, Secretaría Distrital de Movilidad, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 90),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

Dist_ZAT <- viajes_dist[viajes_dist$modo_principal_comparado == "Auto",c("zat_hg", "modo_principal_comparado", "Distancia_incl_caminata", "fexp_vj")]%>%
  left_join(ZAT[,c("ZAT")], by = c("zat_hg" = "ZAT")) %>%
  group_by(zat_hg) %>%
  summarize(dist_auto = sum(fexp_vj*Distancia_incl_caminata/1000))

ZAT <- ZAT %>% left_join(Dist_ZAT, by = c("ZAT" = "zat_hg"))

ZAT$dist_auto[is.na(ZAT$dist_auto)] <- 0

UTAM_indic <- ZAT %>% 
  st_drop_geometry() %>%
  group_by(UTAM) %>% 
  summarise(dist_auto = sum(dist_auto), 
           auto_intensity = sum(dist_auto)/sum(pop))

UTAM <- UTAM %>% left_join(UTAM_indic, by = "UTAM")

bks_car <- round(as.numeric(quantile(UTAM$auto_intensity, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- c(bks_car,max(UTAM$auto_intensity))

UTAM$binned_auto_intensity <- cut(UTAM$auto_intensity,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = binned_auto_intensity))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Límite de municipio fuera de Bogotá DC \nLímite de localidad dentro de Bogotá DC"), color = "grey")+
  #geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 4.5, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(6,"RdPu")) +
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(title = "", fill = "Distancia cotidiana promedio per \ncápita en auto particular (km)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, Secretaría Distrital de Movilidad, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 90),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

Dist_ZAT <- viajes_dist[viajes_dist$modo_principal_comparado == "Transporte publico",c("zat_hg", "modo_principal_comparado", "Distancia_incl_caminata", "fexp_vj")]%>%
  left_join(ZAT[,c("ZAT")], by = c("zat_hg" = "ZAT")) %>%
  group_by(zat_hg) %>%
  summarize(dist_transit = sum(fexp_vj*Distancia_incl_caminata/1000))

ZAT <- ZAT %>% left_join(Dist_ZAT, by = c("ZAT" = "zat_hg"))

ZAT$dist_transit[is.na(ZAT$dist_transit)] <- 0

UTAM_indic <- ZAT %>% 
  st_drop_geometry() %>%
  group_by(UTAM) %>% 
  summarise(dist_transit = sum(dist_transit), 
           transit_intensity = sum(dist_transit)/sum(pop))

UTAM <- UTAM %>% left_join(UTAM_indic, by = "UTAM")

bks_transit <- round(as.numeric(quantile(UTAM$transit_intensity, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),3)
bks <- c(bks_transit,max(UTAM$transit_intensity))

UTAM$binned_transit_intensity <- cut(UTAM$transit_intensity,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = binned_transit_intensity))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Límite de municipio fuera de Bogotá DC \nLímite de localidad dentro de Bogotá DC"), color = "grey")+
  #geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 4.5, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(6,"Greens")) +
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(title = "", fill = "Distancia cotidiana promedio per \ncápita en transporte público (km)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, Secretaría Distrital de Movilidad, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 90),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

Dist_ZAT <- viajes_dist[,c("zat_hg", "modo_principal_comparado", "CO2-eq", "fexp_vj")]%>%
  left_join(ZAT[,c("ZAT")], by = c("zat_hg" = "ZAT")) %>%
  group_by(zat_hg) %>%
  summarize(co2 = sum(fexp_vj*`CO2-eq`/1000000))

ZAT <- ZAT %>% left_join(Dist_ZAT, by = c("ZAT" = "zat_hg"))

ZAT$co2[is.na(ZAT$co2)] <- 0

UTAM_indic <- ZAT %>% 
  st_drop_geometry() %>%
  group_by(UTAM) %>% 
  summarise(co2 = sum(co2), 
           co2_capita = 1000*sum(co2)/sum(pop))

UTAM <- UTAM %>% left_join(UTAM_indic, by = "UTAM")

bks_co2 <- round(as.numeric(quantile(UTAM$co2_capita, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),5)
bks <- c(bks_co2,max(UTAM$co2_capita))

UTAM$binned_co2_capita <- cut(UTAM$co2_capita,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = binned_co2_capita))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Límite de municipio fuera de Bogotá DC \nLímite de localidad dentro de Bogotá DC"), color = "grey")+
  #geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 4.5, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(6,"Oranges")) +
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(title = "", fill = "Emisiones cotidianas promedio \nper cápita (kg CO2-eq)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, Secretaría Distrital de Movilidad, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 90),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

Cartography in English

#Colores
bks_total <- round(as.numeric(quantile(UTAM$dist_capita, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),2)
bks <- c(bks_total,max(UTAM$dist_capita))

UTAM$binned_dist_capita <- cut(UTAM$dist_capita,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = binned_dist_capita))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 4.5, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(5,"Purples")) +
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(title = "", fill = "Average daily distance \nper capita (km)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Bogotá Mobility Department (SDM), 2023 \nDiscretization in quintiles", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

#Colores
bks_walk <- round(as.numeric(quantile(UTAM$walk_intensity, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),2)
bks <- c(bks_walk,max(UTAM$walk_intensity))

UTAM$binned_walk_intensity <- cut(UTAM$walk_intensity,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = binned_walk_intensity))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 2.5, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(5,"Reds")) +
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(title = "", fill = "Average daily walking distance \nper capita (km)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Bogotá Mobility Department (SDM), 2023 \nDiscretization in quintiles", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

bks_car <- round(as.numeric(quantile(UTAM$auto_intensity, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),2)
bks <- c(bks_car,max(UTAM$auto_intensity))

UTAM$binned_auto_intensity <- cut(UTAM$auto_intensity,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = binned_auto_intensity))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 2.5, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(5,"RdPu")) +
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(title = "", fill = "Average daily private car distance \nper capita (km)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Bogotá Mobility Department (SDM), 2023 \nDiscretization in quintiles", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

bks_transit <- round(as.numeric(quantile(UTAM$transit_intensity, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),3)
bks <- c(bks_transit,max(UTAM$transit_intensity))

UTAM$binned_transit_intensity <- cut(UTAM$transit_intensity,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = binned_transit_intensity))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 2.5, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(5,"Greens")) +
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
  labs(title = "", fill = "Average daily public transport distance \nper capita (km)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Bogotá Mobility Department (SDM), 2023 \nDiscretization in quintiles", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

bks_co2 <- round(as.numeric(quantile(UTAM$co2_capita, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),5)
bks <- c(bks_co2,max(UTAM$co2_capita))

UTAM$binned_co2_capita <- cut(UTAM$co2_capita,bks)

ggplot()+
  theme_bw()+
  geom_sf(data = UTAM, color = NA, aes(fill = binned_co2_capita))+
  geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
  #geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 2.5, check_overlap = FALSE) +
  scale_fill_manual(values = brewer.pal(5,"Oranges")) +
  coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
 labs(title = "", fill = "Average daily GHG emissions \nper capita (kg CO2-eq)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Bogotá Mobility Department (SDM), 2023 \nDiscretization in quintiles", linetype = "") +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
        legend.title = element_text(size = 13),
        legend.text = element_text(size=13),
        plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
        legend.text.align = 1)+
  guides(linetype = guide_legend(label.hjust = 0, order = 1))+
  labs(x = "", y = "") +
  annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
  annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))

Correlaciones

setwd("C:/Users/hugot/Documents/Uniandes/LATS")

UTAM_ICS <- read.xlsx("UTAM_ICS.xlsx")

UTAM <- UTAM %>% left_join(UTAM_ICS[, c("UTAM", "propICS12")], by = "UTAM")

correlaciones <- rbind(data.frame(cor.test(UTAM$propICS12, UTAM$auto_intensity, method = "pearson")[c(4,3)]),
              data.frame(cor.test(UTAM$walk_intensity, UTAM$auto_intensity, method = "pearson")[c(4,3)]),
              data.frame(cor.test(UTAM$walk_intensity, UTAM$transit_intensity, method = "pearson")[c(4,3)]),
              data.frame(cor.test(UTAM$propICS12, UTAM$walk_intensity, method = "pearson")[c(4,3)]),
              data.frame(cor.test(UTAM$propICS12, UTAM$transit_intensity, method = "pearson")[c(4,3)]),
              data.frame(cor.test(UTAM$co2_capita, UTAM$auto_intensity, method = "pearson")[c(4,3)]),
              data.frame(cor.test(UTAM$propICS12, UTAM$co2_capita, method = "pearson")[c(4,3)]))

rownames(correlaciones) <- c("auto_ICS",
                             "walk_auto",
                             "walk_transit",
                             "ICS_walk",
                             "ICS_transit",
                             "GES_auto",
                             "ICS_GES")

kable(correlaciones, caption = "Correlaciones", format.args = list(big.mark = ","))

Correlaciones
	estimate	p.value
auto_ICS	-0.4873969	0.0000000
walk_auto	-0.4482729	0.0000001
walk_transit	0.4210982	0.0000005
ICS_walk	0.4398057	0.0000001
ICS_transit	0.5627718	0.0000000
GES_auto	0.8709694	0.0000000
ICS_GES	-0.2907798	0.0007185

correlaciones <- cbind(rownames(correlaciones), correlaciones)

write.xlsx(correlaciones, "correlaciones.xlsx", overwrite = TRUE)

Palma ratio

UTAM_top_10 <- UTAM[order(UTAM$propICS12, decreasing = FALSE),]
select <-as.integer(0.1*nrow(UTAM_top_10))
UTAM_top_10 <- UTAM_top_10[1:select,]

UTAM_bottom_40 <- UTAM[order(UTAM$propICS12, decreasing = TRUE),]
select <-as.integer(0.4*nrow(UTAM_bottom_40))
UTAM_bottom_40 <- UTAM_bottom_40[1:select,]

palma <- rbind(data.frame(top_10 = weighted.mean(x = UTAM_top_10$dist_capita, w = UTAM_top_10$pop),
                   bottom_40 = weighted.mean(x = UTAM_bottom_40$dist_capita, w = UTAM_bottom_40$pop)),
               data.frame(top_10 = weighted.mean(x = UTAM_top_10$auto_intensity, w = UTAM_top_10$pop),
                   bottom_40 = weighted.mean(x = UTAM_bottom_40$auto_intensity, w = UTAM_bottom_40$pop)),
               data.frame(top_10 = weighted.mean(x = UTAM_top_10$walk_intensity, w = UTAM_top_10$pop),
                   bottom_40 = weighted.mean(x = UTAM_bottom_40$walk_intensity, w = UTAM_bottom_40$pop)),
               data.frame(top_10 = weighted.mean(x = UTAM_top_10$transit_intensity, w = UTAM_top_10$pop),
                   bottom_40 = weighted.mean(x = UTAM_bottom_40$transit_intensity, w = UTAM_bottom_40$pop)),
               data.frame(top_10 = weighted.mean(x = UTAM_top_10$co2_capita, w = UTAM_top_10$pop),
                   bottom_40 = weighted.mean(x = UTAM_bottom_40$co2_capita, w = UTAM_bottom_40$pop)))


palma$palma <- palma$top_10/palma$bottom_40

rownames(palma) <- c("dist_capita",
                     "auto_intensity",
                     "walk_intensity",
                     "transit_intensity",
                     "co2_capita")

kable(palma, caption = "Palma ratios", format.args = list(big.mark = ","))

Palma ratios
	top_10	bottom_40	palma
dist_capita	10.4797194	14.1213597	0.7421183
auto_intensity	3.6798541	1.0520899	3.4976612
walk_intensity	0.7930449	1.4220469	0.5576784
transit_intensity	3.5481285	8.6480176	0.4102823
co2_capita	0.8639931	0.6056498	1.4265556

palma <- cbind(rownames(palma), palma)

write.xlsx(palma, "palma.xlsx", overwrite = TRUE)

Valores promediados de los indicadores

promedios <- data.frame(dist_capita = sum(UTAM$dist)/sum(UTAM$pop),
                        walk_intensity = sum(UTAM$dist_walk)/sum(UTAM$pop),
                        auto_intensity = sum(UTAM$dist_auto)/sum(UTAM$pop),
                        transit_intensity = sum(UTAM$dist_transit)/sum(UTAM$pop),
                        ges_capita = 1000*sum(UTAM$co2)/sum(UTAM$pop))

kable(promedios, caption = "Valores promediados de los indicadores", format.args = list(big.mark = ","))

Valores promediados de los indicadores
dist_capita	walk_intensity	auto_intensity	transit_intensity	ges_capita
12.71442	1.161807	1.72813	7.015431	0.6600684

Distancias Bogotá EODH 2023

Hugo Thomas

2025-03-06